Нашей задачей было построение модели для вычисления возраста опопссумов по их морфологическим данным. Исходя из литературного поиска, мы предполагаем, что используемый набор данных был получен группой австралийских ученых в октябре-ноябре 1993 года [1, 2]. Для данного исследований были отловлены и измерены взрослые особи Trichosurus caninus, горного короткоухого поссума.
Статистическая обработка проводилась при помощи R (версия).
Загрузка необходимых пакетов производилась с помощью функции, обнаруживающей уже установленные пакеты и определяющей необходимость их загрузки:
package_installer <- function(package){
if (!require(package, character.only=T, quietly=T)) {
install.packages(package)
library(package, character.only=T)
}
}
lapply(c(library(ggplot2),
library(dplyr),
library(ggcorrplot),
library(psych),
library(stats),
library(car),
library(DAAG),
library(cowplot),
library(bslib),
library(gridExtra)), package_installer)
В R данные являются составной частью пакета DAAG. Данные содержат следующие переменные (site — одна из семи локаций отлова поссумов: Cambarville, Bellbird, Whian Whian, Byrangery, Conondale, Allyn River и Bulburin соответственно; Pop – одна из трех популяций поссумов: Victoria, New South Wales и Queensland; sex — факторная переменная с уровнями f - самки, m - самцы; age — возраст; hdlngth — длина головы; skull — ширина черепа; totlngth — общая длина; taill — длина хвоста; footlgth — длина стопы; earconch — длина ушной раковины; eye — расстояние от медиального угла глазной щели до латерального угла глазной щели правого глаза; chest — обхват груди (в см); belly — обхват живота (в см).
Были произведены перекодировка факторных переменных и удаление переменной “Pop” (по требованию заказчика).
possum <- DAAG::possum[-3]
possum$sex <- factor(possum$sex, levels = c("m", "f"), labels = c("male",
"female"))
possum$site <- factor(possum$site, levels = 1:7, labels = c("Cambarville",
"Bellbird",
"Whian Whian",
"Byrangery",
"Conondale",
"Allyn River",
"Bulburin"))
Краткая сводка по данным:
str(possum)
## 'data.frame': 104 obs. of 13 variables:
## $ case : num 1 2 3 4 5 6 7 8 9 10 ...
## $ site : Factor w/ 7 levels "Cambarville",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ sex : Factor w/ 2 levels "male","female": 1 2 2 2 2 2 1 2 2 2 ...
## $ age : num 8 6 6 6 2 1 2 6 9 6 ...
## $ hdlngth : num 94.1 92.5 94 93.2 91.5 93.1 95.3 94.8 93.4 91.8 ...
## $ skullw : num 60.4 57.6 60 57.1 56.3 54.8 58.2 57.6 56.3 58 ...
## $ totlngth: num 89 91.5 95.5 92 85.5 90.5 89.5 91 91.5 89.5 ...
## $ taill : num 36 36.5 39 38 36 35.5 36 37 37 37.5 ...
## $ footlgth: num 74.5 72.5 75.4 76.1 71 73.2 71.5 72.7 72.4 70.9 ...
## $ earconch: num 54.5 51.2 51.9 52.2 53.2 53.6 52 53.9 52.9 53.4 ...
## $ eye : num 15.2 16 15.5 15.2 15.1 14.2 14.2 14.5 15.5 14.4 ...
## $ chest : num 28 28.5 30 28 28.5 30 30 29 28 27.5 ...
## $ belly : num 36 33 34 34 33 32 34.5 34 33 32 ...
summary(possum)
## case site sex age hdlngth
## Min. : 1.00 Cambarville:33 male :61 Min. :1.000 Min. : 82.50
## 1st Qu.: 26.75 Bellbird :13 female:43 1st Qu.:2.250 1st Qu.: 90.67
## Median : 52.50 Whian Whian: 7 Median :3.000 Median : 92.80
## Mean : 52.50 Byrangery : 7 Mean :3.833 Mean : 92.60
## 3rd Qu.: 78.25 Conondale :13 3rd Qu.:5.000 3rd Qu.: 94.72
## Max. :104.00 Allyn River:13 Max. :9.000 Max. :103.10
## Bulburin :18 NA's :2
## skullw totlngth taill footlgth
## Min. :50.00 Min. :75.00 Min. :32.00 Min. :60.30
## 1st Qu.:54.98 1st Qu.:84.00 1st Qu.:35.88 1st Qu.:64.60
## Median :56.35 Median :88.00 Median :37.00 Median :68.00
## Mean :56.88 Mean :87.09 Mean :37.01 Mean :68.46
## 3rd Qu.:58.10 3rd Qu.:90.00 3rd Qu.:38.00 3rd Qu.:72.50
## Max. :68.60 Max. :96.50 Max. :43.00 Max. :77.90
## NA's :1
## earconch eye chest belly
## Min. :40.30 Min. :12.80 Min. :22.0 Min. :25.00
## 1st Qu.:44.80 1st Qu.:14.40 1st Qu.:25.5 1st Qu.:31.00
## Median :46.80 Median :14.90 Median :27.0 Median :32.50
## Mean :48.13 Mean :15.05 Mean :27.0 Mean :32.59
## 3rd Qu.:52.00 3rd Qu.:15.72 3rd Qu.:28.0 3rd Qu.:34.12
## Max. :56.20 Max. :17.80 Max. :32.0 Max. :40.00
##
Для удобства дальнейшей работы с данными числовые переменные были выделены в отдельный сабсет.
possum_numerical <- possum[4:13]
Подсчет числа NA:
NA_count <- sum(is.na(possum))
NA - 3 наблюдения из 104 - менее 3%, они могут быть удалены:
possum <- na.omit(na.omit(possum))
possum_numerical <- na.omit(na.omit(possum_numerical))
Проверка на нормальность распределения данных необходима, чтобы оценить круг методов, которые могут применять при дальнейшем анализе данных. При несоблюдении условия нормальности как по формальным критериям (тесты Шапиро-Уилка и Колмогорова-Смирнова), так и по визуальным характеристикам распределения (график плотности и квантиль-квантильный график) используются непараметрические статистические критерии. Они уступают по мощности параметрическим критериям, однако последние применимы только при нормальности распределения.
density <- apply(possum_numerical, 2, function(a) ggplot(possum_numerical, aes(x=a))+
geom_density())
График плотности для интересующей переменной вызывается при помощи команды
density$<variable>, где<variable>– имя переменной. Пример -density$footlgth. Дополнительную кастомизацию можно провести, работая, как с обычной переменной, хранящей график ggplot.
График плотности для возраста поссумов. Распределение скошено влево, не является нормальным.
График плотности для длины головы поссумов. Распределение близко к нормальному. Имеется небольшой субпик в левой части графика.
График плотности для ширины черепа поссумов. Распределение скошено влево, в правой части графика имеется небольшой субпик.
График плотности для длины тела поссумов. Распределение близко к нормальному, с несколько скошенной в правую сторону вершиной.
График плотности для длины хвоста поссумов. Распределение близко к нормальному, несколько бимодально.
График плотности для длины стоп поссумов. Распределение бимодально.
График плотности для длины ушной раковины поссумов. Распределение бимодально.
График плотности для ширины глазной щели поссумов. Распределение близко к нормальному.
График плотности для обхвата груди поссумов. Распределение близко к нормальному, вершина несколько скошено.
График плотности для обхвата живота поссумов. Распределение близко к нормальному.
Квантиль-квантильный график для интересующей переменной вызывается при помощи команды
qqplots$<variable>, где<variable>– имя переменной. Пример -qqplots$earconch. Дополнительную кастомизацию можно провести, работая, как с обычной переменной, хранящей график ggplot.
Квантиль-квантильный график для возраста поссумов. Квантили указывают на распределение, отличающиеся от нормального.
Квантиль-квантильный график для длины головы поссумов. Квантили распределены нормально.
Квантиль-квантильный график для ширины головы поссумов. Распределение квантилей отличается от нормального.
Квантиль-квантильный график для длины тела поссумов. Распределение квантилей близко к нормальному.
Квантиль-квантильный график для длины хвоста поссумов. Распределние квантилей визуально несколько отличается от нормального.
Квантиль-квантильный график для длины стоп поссумов. Распределение квантилей отличается от нормального.
Квантиль-квантильный график для длины ушной раковины поссумов. Распределение квантилей отличается от нормального.
Квантиль-квантильный график для ширины глазничной щели поссумов. Распределение квантилей близко к нормальному.
Квантиль-квантильный график для обхвата груди поссумов. Распределение квантилей близко к нормальному.
Квантиль-квантильный график для обхвата живота поссумов. Распределение квантилей близко к нормальному.
Следующим шагом проверки на нормальность стало использование теста Шапиро-Уилка для указанных данных.
normality <- sapply(possum_numerical, function(a) shapiro.test(a))
normality
## age hdlngth
## statistic 0.9374365 0.9813663
## p.value 0.0001246421 0.1648769
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "a" "a"
## skullw totlngth
## statistic 0.9381717 0.9848761
## p.value 0.0001380455 0.3045355
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "a" "a"
## taill footlgth
## statistic 0.9830628 0.9518154
## p.value 0.2228095 0.001023365
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "a" "a"
## earconch eye
## statistic 0.9171524 0.9777459
## p.value 9.083728e-06 0.08554447
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "a" "a"
## chest belly
## statistic 0.9850221 0.9888213
## p.value 0.3121142 0.5633678
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "a" "a"
По результатам анализа графиков и теста Шапиро-Уилка можно сделать вывод, что длина головы, общая длина, длина хвоста, ширина глазной щели, обхват груди и обхват живота распределены нормально (p-value теста Шапиро-Уилка больше 0.05, графики плотности и квантилей близки к нормальным), а возраст, ширина черепа, длина стоп и длина ушной раковины - ненормально.
Боксплоты по количественным переменным:
boxplot(possum$age, possum$hdlngth, possum$skullw, possum$totlngth,
possum$taill, possum$footlgth, possum$earconch, possum$eye,
possum$chest, possum$belly,
names = c('age',
'hdlngth',
'skullw',
'totlngth',
'taill',
'footlgth',
'earconch',
'eye',
'chest',
'belly'))
В переменных hdlngth, skullw, taill, eye, chest, belly отмечены выбросы в количестве, не допускающем их удаление. Далее можно проанализировать их и сделать предположения о биологическом смысле некоторых из них.
Чтобы при подборе модели избежать возможных проблем с коллинеарностью, был проведен корреляционный анализ данных. Была составлена матрица корреляций всех колчиественных переменных. Корреляции с p-value меньше 0.05 не считались статистически значимыми и помечались на графике крестом. Поправка на множественные сравнения не проводилась, так как результаты анализа влияли лишь на ужесточение требование к модели, а не на непосредственные выдвижение и проверку гипотез.
pleiades_mat <- corr.test(possum_numerical)
temperature_plot <- ggcorrplot(pleiades_mat$r, type = "lower", outline.col = "white", lab=TRUE, method = "circle", p.mat = pleiades_mat$p)
temperature_plot
Мы наблюдали сильную положительную связь между длиной и шириной черепа, а также между длиной стоп и ушных раковин. Мы обнаружили положительную связь средней силы между обхватом поясницы и обхватом груди, обхватом поясницы и общей длиной тела, обхватом поясницы и длиной головы, обхватом груди и длиной тела, обхватом груди и шириной черепа, обхватом груди и высотой черепа, общей длиной тела и высотой черепа. Мы обнаружили умеренную отрицательную связь между длиной хвоста и длиной ушной раковины. Также мы наблюдали умеренную положительную связь между обхватом поясницы и длиной стоп, обхватом пояницы и шириной черепа, обхватом поясницы и возрастом, обхватом груди и длиной стоп, обхватом груди и возрастом, шириной глазной щели и шириной черепа, шириной глазной щели и длиной головы, длиной стоп и общей длиной тела, длиной стоп и длиной головы, возрастом и длиной головы.
Большинство указанных корреляций было ожидаемо по биологическим причинам (пропорцианольность строения тела и его отдельных частей). Наше внимание привлекла корреляция между длиной ушной раковины и длиной стоп.
Для дальнейшего анализа были построены попарные диаграммы рассеивания. На тех из них, которые содержали переменные footlhth и earconch обнаруживались два кластера. Цвет точек в первом наборе диаграмм обозначал пол, во втором - локацию поссума.
deviant_apply_sex <- apply(possum_numerical, 2, function(b) apply(possum_numerical, 2, function(a) ggplot(possum)+
geom_point(aes(x=a, y=b, col=sex))))
deviant_apply_site <- apply(possum_numerical, 2, function(b) apply(possum_numerical, 2, function(a) ggplot(possum)+
geom_point(aes(x=a, y=b, col=site))))
Диаграмма рассеивания для интересующих переменных вызывается при помощи команды
density$<variable_y>$<variable_x>, где<variable_y>– имя переменной, которая будет располагаться по оси ординат, а<variable_x>– имя переменной, которая будет располагаться по оси абсцисс. Пример -density$footlgth. Дополнительную кастомизацию можно провести, работая, как с обычной переменной, хранящей график ggplot.
Для экономии места приведена будет наиболее показательная диаграмма рассеивания между переменными footlhth и earconch, с указанием локации измерения поссума.
Заметно, что поссумы с более длинными стопами и ушными раковинами принадлежат к популяции, обитающей в первом и втором локусе, а с более маленькими – в остальных пяти. Таким образом, можно предположить, что мы имеем дело с двумя подвидами поссумов, отличающихся по морфометрическим характеристикам и ареалам обитания.
Для более подного разведочного анализа было решено оценить распределения переменных при разбиении их на подгруппы по полу и локации.
possum$locus <- factor(ifelse(possum$site %in% c("Cambarville", "Bellbird"), "First_locus", "Second_locus"))
possum_numerical$locus <- possum$locus
Добавленная переменная соответствует гипотетическим подвидам, обитающим в одном из двух локусов.
Распределения асимметричны, скошены влево.
Распределения в целом выглядят как нормальные
Заметная ассимметрия (скошены влево)
Распределение немного скошено вправо в случае самок, имеется два пика
у самцов
Распределения не очень симметричны, в целом выглядят как нормальные
Распределения сильно отличаются от нормальных: имеют по два пика,
форму, сильно отличную от нормального распределения
Распределения сильно отличаются от нормальных: имеют по два пика,
форму, сильно отличную от нормального распределения
Распределения не очень симметричны, в целом выглядят как нормальные
У распределений немного искажена форма, в целом выглядят приближенно
к нормальным
Распределения близки к нормальным. У поссумов женского пола имеется малый пик в правой части графика. Было высказано предположение о том, что это беременные самки.
Распределения скошены влево.
Распределения в целом выглядят как нормальные
Распределения близки к нормальным.
Распределение у поссумов из второго локуса близко к нормальному, у поссумов первого локуса несколько скошено вправо.
Распределения нормальны, заметно различие между двумя популяциями.
Распределения нормальны, заметно различие между двумя популяциями.
Распределения нормальны, заметно различие между двумя популяциями.
Распределение в первому локусе близко к нормальному, во втором имеет утолщенный хвост справа.
Распределения асимметричны.
Распределение в первом локусе близко к нормальному, во втором имеет утолщенные хвосты.
При статистическом анализе морфологических данных мы обнаружили выраженное бимодальное распределение по некоторым показателям (размер ушной раковины, стопы и длина хвоста). Для объяснения данной гетерогенности популяции был проведен литературный поиск, который выявил, что спустя 9 лет после сбора исследуемых данных о Trichosurus caninus, в 2002 году, было предложено разделить этот вид на северный (сохранил название Trichosurus caninus) и южный (T. cunninghami) подвиды, вследствие существенных морфологических отличий (при этом генетических оснований для разделения на отдельные виды оказалось недостаточно)[2]. То есть, данные о горном короткоухом поссуме, собранные до 2002 года, объединяют в себе сведения о двух морфологически отличных подвидах.
На рисунке 1 продемонстрированы морфологические различия этих двух видов [2]. Согласно данным австралийских коллег, у южных популяций: – более крупная ушная раковина, – значительно (P < 0,001) более длинная стопа и – значительно (P < 0,001) более короткий хвост, чем у особей из северных популяций. Животных можно четко различить с помощью индекса, основанного на этих трех морфологических показателях. Эти данные полностью согласуются с результатами нашего исследования.
Кроме того разброс значений по некоторым показателям может быть вызван примесью особей кистехвостых поссумов (Trichosurus vulpecula), поскольку они имеют схожую морфологию и пересечение ареала обитания с Trichosurus caninus.
Рисунок 1. Ключевые различия в морфологии T. cuninghami и T. caninus (по D. B. Lindenmayer et al.)
Нами также была обнаружена незначительная бимодальность в распределении величины обхвата живота: 4 наблюдения образуют вторую моду в области больших значений. Учитывая, что эти данные описывают женские особи (среди мужских особей второй моды у данного распределения не обнаружено), мы предположили, что они относятся к беременным / вынашивающим самкам. Отсутствие больших значений по этому показателю среди юных самок также согласуется с данными о возрасте половой зрелости, который наступает на третьем году жизни.
Анализ литературных данных о беременности горных поссумов показал, что внутриутробная беременность Trichosurus caninus длится чуть больше двух недель и приводит к рождению недоношенного плода очень малого размера, что вероятно, если и приводит, то к незначительному увеличению живота.
После рождения детёныш перебирается в сумку матери, где безвылазно проводит следующие 175–200 дней. Поскольку детёныш на протяжении этого срока является неотъёмным от матери, замеры живота кормящих самок либо а) проводятся с учётом младенца, либо б) отбрасываются (если дизайн эксперимента позволяет отобрать особи среди некормящих самок).
На следующем этапе вскармливания младенец больше времени проводит на спине у матери, но при этом продолжает посещать сумку до полного отъёма от груди в возрасте от 240 до 270 дней. То есть животы самок и на этом сроке выкармливания могут быть существенно увеличены за счет детеныша
Участвовали ли кормящие самки в данном исследовании не уточняется, однако, согласно данным другого исследования (проводившегося в тот же год в тех же регионах), доля вынашивающих самок в популяции составляла около 88% в летний период и 60% к концу осени (Таблица 1, [3]).
Учитывая такой высокий процент среди женской популяции, разумно предположить, что кормящие самки не были отброшены при проведении замеров. Однако вторая мода в распределении величины живота образована намного меньшей долей замеров, что не согласуется с ожидаемым высоким процентом вынашивающих самок, и может быть вызвана другими обстоятельствами: наличием нескольких детенышей в либо вариабельностью их размеров.
Таблица 1. Доля самок с детенышами (по Viggers K.L. с изменениями)
| Апрель | Июнь | Август | Ноябрь |
|---|---|---|---|
| 0.25(2/8) | 0.88 (7/8) | 0.88 (7/8) | 0.60 (6/10) |
Средняя продолжительности жизни у обоих подвидов Trichosurus caninus составляет 8-12 лет, при этом встречаются особи возрастом 17 лет [1]. Поссумы этого вида считаются долгожителями среди сородичей (к примеру, упомянутый выше кистехвостый Trichosurus vulpecula живёт в среднем 7-10 лет). Однако в представленных данных медианный возраст особи составил всего 3 года, а максимальный – 9 лет.
Поскольку сбор данных производился путём отлова животных в естественной среде обитания (а не в контролируемом непрерывном наблюдении), их возраст доподлинно неизвестен и был оценен косвенным способом – по степени стираемости зубов. В отличие от других измерений, полученных прямым способом (замеры частей тела).
При таком методе оценки может возникать дополнительная погрешность, вызванная травмами челюсти или генетическими дефектами эмали или пищеварения и др. Это может негативно сказаться на построении предсказательной модели.
Модели по одному предиктору:
mod <- lm(age~hdlngth, data=possum)
summary(mod)
##
## Call:
## lm(formula = age ~ hdlngth, data = possum)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.6815 -1.3815 -0.2424 1.1479 5.0761
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -12.80894 4.79274 -2.673 0.008804 **
## hdlngth 0.17934 0.05165 3.472 0.000766 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.817 on 99 degrees of freedom
## Multiple R-squared: 0.1086, Adjusted R-squared: 0.09957
## F-statistic: 12.06 on 1 and 99 DF, p-value: 0.0007661
mod <- lm(age~skullw, data=possum)
summary(mod)
##
## Call:
## lm(formula = age ~ skullw, data = possum)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.715 -1.441 -0.300 1.347 5.295
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.21855 3.39152 -1.834 0.0697 .
## skullw 0.17627 0.05945 2.965 0.0038 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.845 on 99 degrees of freedom
## Multiple R-squared: 0.08155, Adjusted R-squared: 0.07227
## F-statistic: 8.79 on 1 and 99 DF, p-value: 0.003795
mod <- lm(age~totlngth, data=possum)
summary(mod)
##
## Call:
## lm(formula = age ~ totlngth, data = possum)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2173 -1.4215 -0.1766 1.1705 4.9051
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.86308 3.86020 -1.778 0.07849 .
## totlngth 0.12244 0.04418 2.771 0.00667 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.854 on 99 degrees of freedom
## Multiple R-squared: 0.07198, Adjusted R-squared: 0.06261
## F-statistic: 7.679 on 1 and 99 DF, p-value: 0.006673
mod <- lm(age~taill, data=possum)
summary(mod)
##
## Call:
## lm(formula = age ~ taill, data = possum)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0495 -1.2322 -0.5825 1.4175 5.1840
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.50415 3.59572 -0.140 0.889
## taill 0.11676 0.09692 1.205 0.231
##
## Residual standard error: 1.911 on 99 degrees of freedom
## Multiple R-squared: 0.01445, Adjusted R-squared: 0.004494
## F-statistic: 1.451 on 1 and 99 DF, p-value: 0.2312
mod <- lm(age~footlgth, data=possum)
summary(mod)
##
## Call:
## lm(formula = age ~ footlgth, data = possum)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0847 -1.5317 -0.5043 1.3643 4.9591
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.07640 2.96524 0.026 0.979
## footlgth 0.05476 0.04326 1.266 0.209
##
## Residual standard error: 1.909 on 99 degrees of freedom
## Multiple R-squared: 0.01592, Adjusted R-squared: 0.005984
## F-statistic: 1.602 on 1 and 99 DF, p-value: 0.2086
mod <- lm(age~earconch, data=possum)
summary(mod)
##
## Call:
## lm(formula = age ~ earconch, data = possum)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.9926 -1.6770 -0.6708 1.2980 5.0793
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.31802 2.28480 1.015 0.313
## earconch 0.03124 0.04730 0.660 0.510
##
## Residual standard error: 1.921 on 99 degrees of freedom
## Multiple R-squared: 0.004387, Adjusted R-squared: -0.00567
## F-statistic: 0.4362 on 1 and 99 DF, p-value: 0.5105
mod <- lm(age~eye, data=possum)
summary(mod)
##
## Call:
## lm(formula = age ~ eye, data = possum)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2201 -1.4650 -0.3812 1.1994 5.2413
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.4912 2.6684 -0.934 0.3528
## eye 0.4195 0.1769 2.372 0.0196 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.872 on 99 degrees of freedom
## Multiple R-squared: 0.05376, Adjusted R-squared: 0.0442
## F-statistic: 5.624 on 1 and 99 DF, p-value: 0.01965
mod <- lm(age~chest, data=possum)
summary(mod)
##
## Call:
## lm(formula = age ~ chest, data = possum)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.7539 -1.2776 -0.1663 0.8811 4.8811
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.77202 2.43571 -1.959 0.052903 .
## chest 0.31753 0.08975 3.538 0.000616 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.814 on 99 degrees of freedom
## Multiple R-squared: 0.1122, Adjusted R-squared: 0.1033
## F-statistic: 12.52 on 1 and 99 DF, p-value: 0.0006156
mod <- lm(age~belly, data=possum)
summary(mod)
##
## Call:
## lm(formula = age ~ belly, data = possum)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.913 -1.407 -0.420 1.340 5.087
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.44667 2.15541 -2.063 0.04173 *
## belly 0.25333 0.06581 3.849 0.00021 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.795 on 99 degrees of freedom
## Multiple R-squared: 0.1302, Adjusted R-squared: 0.1214
## F-statistic: 14.82 on 1 and 99 DF, p-value: 0.00021
В качестве единственных предикторов значимо предсказывают возраст переменные chest (***), belly (***), hdlngth (***), skullw (**), totlngth (**), eye (*).
Предсказательная сила простых линейных моделей оказалсь неудовлетворительной, поэтому была предпринята попытка построить множественную линейную модель для предсказания возраст поссумов.
Стандартизируем модель, чтобы получить информацию о том, какие из возможных предикторов оказывают наиболшее влияние на полную модель.
possum_scale <- as.data.frame(sapply(possum_numerical[-11], scale))
mod_scale <- lm(age ~ ., data = possum_scale)
summary(mod_scale)
##
## Call:
## lm(formula = age ~ ., data = possum_scale)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.64661 -0.62064 -0.09882 0.49549 2.60641
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.804e-16 9.364e-02 0.000 1.000
## hdlngth 8.912e-02 1.648e-01 0.541 0.590
## skullw 3.445e-02 1.437e-01 0.240 0.811
## totlngth 4.928e-02 1.814e-01 0.272 0.787
## taill -3.262e-04 1.479e-01 -0.002 0.998
## footlgth -2.512e-01 1.916e-01 -1.311 0.193
## earconch 2.081e-01 1.808e-01 1.151 0.253
## eye 1.395e-01 1.047e-01 1.333 0.186
## chest 1.595e-01 1.495e-01 1.067 0.289
## belly 2.047e-01 1.273e-01 1.608 0.111
##
## Residual standard error: 0.9411 on 91 degrees of freedom
## Multiple R-squared: 0.1941, Adjusted R-squared: 0.1144
## F-statistic: 2.436 on 9 and 91 DF, p-value: 0.01574
Заметно, что сильнее всего на модель влияют belly, earconch и footlgth.
Удаляем коллинеарные факторы посредством функции вздутия модели.
multiple_mod1 <- lm(age ~., data=possum_numerical[-11])
summary(multiple_mod1)
##
## Call:
## lm(formula = age ~ ., data = possum_numerical[-11])
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.1535 -1.1886 -0.1893 0.9489 4.9917
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.369e+01 6.501e+00 -2.106 0.038 *
## hdlngth 4.851e-02 8.969e-02 0.541 0.590
## skullw 2.127e-02 8.870e-02 0.240 0.811
## totlngth 2.249e-02 8.280e-02 0.272 0.787
## taill -3.169e-04 1.436e-01 -0.002 0.998
## footlgth -1.090e-01 8.316e-02 -1.311 0.193
## earconch 9.813e-02 8.528e-02 1.151 0.253
## eye 2.524e-01 1.894e-01 1.333 0.186
## chest 1.512e-01 1.417e-01 1.067 0.289
## belly 1.437e-01 8.936e-02 1.608 0.111
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.802 on 91 degrees of freedom
## Multiple R-squared: 0.1941, Adjusted R-squared: 0.1144
## F-statistic: 2.436 on 9 and 91 DF, p-value: 0.01574
vif(multiple_mod1)
## hdlngth skullw totlngth taill footlgth earconch eye chest
## 3.0664 2.3316 3.7174 2.4696 4.1471 3.6911 1.2376 2.5238
## belly
## 1.8292
multiple_mod2 <- update(multiple_mod1, .~. - footlgth)
vif(multiple_mod2)
## hdlngth skullw totlngth taill earconch eye chest belly
## 3.0546 2.3233 3.4749 2.4636 1.7309 1.2375 2.4756 1.8234
multiple_mod3 <- update(multiple_mod2, .~. - totlngth)
vif(multiple_mod3)
## hdlngth skullw taill earconch eye chest belly
## 2.6009 2.3231 1.4069 1.3919 1.2333 2.3826 1.8229
multiple_mod4 <- update(multiple_mod3, .~. - hdlngth)
vif(multiple_mod4)
## skullw taill earconch eye chest belly
## 1.8038 1.3720 1.3469 1.1760 2.2964 1.7594
multiple_mod5 <- update(multiple_mod4, .~. - chest)
vif(multiple_mod5)
## skullw taill earconch eye belly
## 1.3433 1.3716 1.2602 1.1610 1.3685
summary(multiple_mod5)
##
## Call:
## lm(formula = age ~ skullw + taill + earconch + eye + belly, data = possum_numerical[-11])
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3800 -1.2481 -0.3102 1.0757 4.8933
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -12.24491 5.72649 -2.138 0.0351 *
## skullw 0.07380 0.06695 1.102 0.2731
## taill 0.01277 0.10646 0.120 0.9048
## earconch 0.03260 0.04956 0.658 0.5123
## eye 0.24818 0.18243 1.360 0.1769
## belly 0.18645 0.07687 2.426 0.0172 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.792 on 95 degrees of freedom
## Multiple R-squared: 0.1679, Adjusted R-squared: 0.1241
## F-statistic: 3.835 on 5 and 95 DF, p-value: 0.003299
После этого шага имеем модель: age = -11.64 + 0.076 * skullw + 0.003 * taill + 0.024 * earconch + 0.259 * eye + 0.183 * belly
Убираем невлияющие предикторы, которые не оказывают значимого влияния на модель.
drop1(multiple_mod5, test = "F")
## Single term deletions
##
## Model:
## age ~ skullw + taill + earconch + eye + belly
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 305.20 123.69
## skullw 1 3.9031 309.10 122.97 1.2149 0.27314
## taill 1 0.0462 305.24 121.70 0.0144 0.90478
## earconch 1 1.3901 306.59 122.15 0.4327 0.51226
## eye 1 5.9459 311.14 123.64 1.8508 0.17691
## belly 1 18.9012 324.10 127.76 5.8834 0.01717 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
multiple_mod6 <- update(multiple_mod5, .~. - taill)
drop1(multiple_mod6, test = "F")
## Single term deletions
##
## Model:
## age ~ skullw + earconch + eye + belly
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 305.24 121.70
## skullw 1 4.0765 309.32 121.05 1.2821 0.26034
## earconch 1 1.4417 306.69 120.18 0.4534 0.50233
## eye 1 5.9882 311.23 121.67 1.8833 0.17316
## belly 1 20.8271 326.07 126.37 6.5502 0.01205 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
multiple_mod7 <- update(multiple_mod6, .~. - earconch)
drop1(multiple_mod7, test = "F")
## Single term deletions
##
## Model:
## age ~ skullw + eye + belly
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 306.69 120.18
## skullw 1 4.2471 310.93 119.57 1.3433 0.249299
## eye 1 5.1815 311.87 119.87 1.6388 0.203540
## belly 1 21.9771 328.66 125.17 6.9510 0.009756 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
multiple_mod8 <- update(multiple_mod7, .~. - skullw)
drop1(multiple_mod8, test = "F")
## Single term deletions
##
## Model:
## age ~ eye + belly
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 310.93 119.57
## eye 1 8.107 319.04 120.17 2.5551 0.113159
## belly 1 36.141 347.07 128.68 11.3909 0.001059 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(multiple_mod8)
##
## Call:
## lm(formula = age ~ eye + belly, data = possum_numerical[-11])
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.9596 -1.2969 -0.4105 1.1155 4.9715
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.76642 2.98114 -2.605 0.01061 *
## eye 0.27726 0.17345 1.598 0.11316
## belly 0.22720 0.06732 3.375 0.00106 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.781 on 98 degrees of freedom
## Multiple R-squared: 0.1523, Adjusted R-squared: 0.135
## F-statistic: 8.803 on 2 and 98 DF, p-value: 0.0003049
На данном шаге имеем модель age = 7.76642 + 0.27726 * eye + 0.22720 * belly
Предсказательная сила модели неудовлетворительная. Дальнейшие шаги проделаны исключительно в учебных целях.
Строим график остатков модели
multiple_mod8_diag <- data.frame(fortify(multiple_mod8), possum_numerical[, c(2:7, 9)])
gg_resid <- ggplot(data = multiple_mod8_diag, aes(x = .fitted, y = .stdresid)) +
geom_point() +
geom_hline(yintercept = 0) +
geom_smooth(method = "lm") +
geom_hline(yintercept = 2, color = "red") +
geom_hline(yintercept = -2, color = "red")
gg_resid
## `geom_smooth()` using formula 'y ~ x'
С моделью не все в порядке. Имеется два наблюдения за пределами 2 стандартных отклонений (влиятельные наблюдения). Ряд измерений в центральной части производят впечатление гетероскедастичности.
Строим график расстояний Кука
cook_model <- ggplot(multiple_mod8_diag, aes(x = 1:nrow(multiple_mod8_diag), y = .cooksd)) +
geom_bar(stat = "identity") +
geom_hline(yintercept = 2, color = "red")
cook_model
График расстояний Кука выглядит приемлемо.
Проверяем нормальность распределения
quantile_model <- qqPlot(multiple_mod8_diag$.stdresid)
quantile_model
## [1] 9 11
Квантильный график выглядит не критично.
res_1 <- gg_resid + aes(x = chest)
res_2 <- gg_resid + aes(x = earconch)
res_3 <- gg_resid + aes(x = taill)
res_4 <- gg_resid + aes(x = totlngth)
res_5 <- gg_resid + aes(x = footlgth)
res_6 <- gg_resid + aes(x = hdlngth)
res_7 <- gg_resid + aes(x = skullw)
grid.arrange(res_1, res_2, res_3, res_4, res_5, res_6, res_7, nrow = 2)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
Несмотря на коллинеарность, возвращаем в модель chest, поскольку наблюдаем явную неучтенную зависимость.
summary(multiple_mod8)
##
## Call:
## lm(formula = age ~ eye + belly, data = possum_numerical[-11])
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.9596 -1.2969 -0.4105 1.1155 4.9715
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.76642 2.98114 -2.605 0.01061 *
## eye 0.27726 0.17345 1.598 0.11316
## belly 0.22720 0.06732 3.375 0.00106 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.781 on 98 degrees of freedom
## Multiple R-squared: 0.1523, Adjusted R-squared: 0.135
## F-statistic: 8.803 on 2 and 98 DF, p-value: 0.0003049
multiple_mod9 <- update(multiple_mod8, .~. + chest)
summary(multiple_mod9)
##
## Call:
## lm(formula = age ~ eye + belly + chest, data = possum_numerical[-11])
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3033 -1.2045 -0.2014 0.9790 4.8330
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.00075 3.27037 -3.058 0.00288 **
## eye 0.28205 0.17211 1.639 0.10450
## belly 0.14694 0.08351 1.760 0.08164 .
## chest 0.17668 0.11036 1.601 0.11263
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.767 on 97 degrees of freedom
## Multiple R-squared: 0.1741, Adjusted R-squared: 0.1486
## F-statistic: 6.817 on 3 and 97 DF, p-value: 0.000323
Несмотря на снижение значимость отдельных предикторов, общий R-квадрат несколько возрастает.
Строим график остатков для новой модели
multiple_mod9_diag <- data.frame(fortify(multiple_mod9), possum_numerical[, c(2:6, 9)])
gg_resid2 <- ggplot(data = multiple_mod9_diag, aes(x = .fitted, y = .stdresid)) +
geom_point() +
geom_hline(yintercept = 0) +
geom_smooth(method = "lm") +
geom_hline(yintercept = 2, color = "red") +
geom_hline(yintercept = -2, color = "red")
gg_resid2
## `geom_smooth()` using formula 'y ~ x'
Распределение остатков стало еще хуже. За границей двух стандартных отклонений стало больше влиятельных наблюдений, а паттерн стал еще более воронкообразным.
Строим график Кука для новой модели
cook_model <- ggplot(multiple_mod9_diag, aes(x = 1:nrow(multiple_mod9_diag), y = .cooksd)) +
geom_bar(stat = "identity") +
geom_hline(yintercept = 2, color = "red")
cook_model
График расстояний Кука выглядит приемлемо.
Строим квантильный график для новой модели.
quantile_model <- qqPlot(multiple_mod9_diag$.stdresid)
quantile_model
## [1] 9 11
Квантильный график выглядит некритично.
Низкая предсказательная сила (скорректированный R-квадрат на уровне 0.15) и неудовлетворительное распределение остатков делают модель, по нашему мнению, непригодной для использования. Нижеследующий раздел (предсказания модели) выполянется в учебных целях, чтобы продемонстрировать, какой подход мы бы применяли, если бы модель была пригодна для дальнейшей работы.
MyData <- data.frame(
belly = seq(min(possum_numerical$belly), max(possum_numerical$belly), length.out = 100),
chest = mean(possum_numerical$chest),
eye = mean(possum_numerical$eye))
Predictions <- predict(multiple_mod9, newdata = MyData, interval = 'confidence')
MyData <- data.frame(MyData, Predictions)
Pl_predict <- ggplot(MyData, aes(x = belly, y = fit)) +
geom_ribbon(alpha = 0.2, aes(ymin = lwr, ymax = upr)) +
geom_line() +
ggtitle("Множественная модель")
Pl_predict
[1] Lindenmayer, D. B., Viggers, K. L., Cunningham, R. B., and Donnelly, C. F. (1995). Morphological variation among populations of the mountain brushtail possum, Trichosurus caninus Ogilby (Phalangeridae, Marsupialia). Australian Journal of Zoology 43, 449–458.
[2] Lindenmayer, D. B., Dubach, J., and Viggers, K. L. (2002). Geographic dimorphism in the mountain brushtail possum (Trichosurus caninus): the case for a new species Australian Journal of Zoology 50, 369–393.
[3] Viggers KL, Lindenmayer DB, Cunningham RB, Donnelly CF. The effects of parasites on a wild population of the Mountain Brushtail Possum (Trichosurus caninus) in south-eastern Australia. Int J Parasitol. 1998 May;28(5):747-55. doi: 10.1016/s0020-7519(98)00022-8. PMID: 9650054.